1 Indroducción



En el presente documento se pretende ir comentando los pasos que se están siguiendo para conseguir predecir la cantidad de agua que llegará a la presa de belesar. Para ello se cuenta con la predicción meteorológica que diariamente se realiza con WRF.

2 Datos



2.1 Histórico DHI



Para la elaboracion del modelo predictivo se cuenta con un Histórico de datos 15 minutales aportados por la presa de Belesar en la que se aportan datos como Aportación, Nivel del embalse, Caudal turbinado, lluvia y algunas otras variables.



DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_DHI_Belesar.RDS'))

DHI



2.1.1 Outliers



Uno de los problemas que nos encontramos a la hora de tratar los datos es la presencia de valores atípicos o “outliers”. modo de ejemplo podemos observar el siguiente gráfico de la aportación… donde se aprecia un valor negativo enorme y erróneo sin ninguna duda.



plot(DHI$`APORTACION (m3/s)`,
     xlab = "Date",
     ylab= "Aportacion [m³/s]")



Investigando sobre la mejor manera de eliminar outliers encontre una función hecha por un fanático. que funciona de la siguiente manera y la verdad que me gusta.



outlierKD <- function(dt, var) {
  var_name <- eval(substitute(var),eval(dt))
  tot <- sum(!is.na(var_name))
  na1 <- sum(is.na(var_name))
  m1 <- mean(var_name, na.rm = T)
  par(mar = rep(2, 4))
  par(mfrow=c(2, 2), oma=c(0,0,3,0))
  boxplot(var_name, main="With outliers")
  hist(var_name, main="With outliers", xlab=NA, ylab=NA)
  outlier <- boxplot.stats(var_name)$out
  mo <- mean(outlier)
  var_name <- ifelse(var_name %in% outlier, NA, var_name)
  boxplot(var_name, main="Without outliers")
  hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
  title("Outlier Check", outer=TRUE)
  na2 <- sum(is.na(var_name))
  message("Outliers identified: ", na2 - na1, " from ", tot, " observations")
  message("Proportion (%) of outliers: ", (na2 - na1) / tot*100)
  message("Mean of the outliers: ", mo)
  m2 <- mean(var_name, na.rm = T)
  message("Mean without removing outliers: ", m1)
  message("Mean if we remove outliers: ", m2)
  
    dt[as.character(substitute(var))] <- invisible(var_name)
    assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
    message("Outliers successfully removed", "\n")
    return(invisible(dt))
   
}



Eliminamos valores negativos y pasamos el filtro para outliers.



DHI$`APORTACION (m3/s)`[which(DHI$`APORTACION (m3/s)`<0)]<- NA 

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 826 from 34552 observations
## Proportion (%) of outliers: 2.39059967585089
## Mean of the outliers: 1075.81016949153
## Mean without removing outliers: 131.482568013429
## Mean if we remove outliers: 108.354577773824
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 317 from 33726 observations
## Proportion (%) of outliers: 0.939927652256419
## Mean of the outliers: 410.450378548896
## Mean without removing outliers: 108.354577773824
## Mean if we remove outliers: 105.488153491574
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 228 from 33409 observations
## Proportion (%) of outliers: 0.682450836600916
## Mean of the outliers: 391.157280701754
## Mean without removing outliers: 105.488153491574
## Mean if we remove outliers: 103.525205991381
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 159 from 33181 observations
## Proportion (%) of outliers: 0.479189897833097
## Mean of the outliers: 380.017798742138
## Mean without removing outliers: 103.525205991381
## Mean if we remove outliers: 102.193901944158
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 168 from 33022 observations
## Proportion (%) of outliers: 0.5087517412634
## Mean of the outliers: 373.650892857143
## Mean without removing outliers: 102.193901944158
## Mean if we remove outliers: 100.805797771961
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 50 from 32854 observations
## Proportion (%) of outliers: 0.152188470201498
## Mean of the outliers: 365.8664
## Mean without removing outliers: 100.805797771961
## Mean if we remove outliers: 100.40179124497
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 4 from 32804 observations
## Proportion (%) of outliers: 0.0121936349225704
## Mean of the outliers: 362.2675
## Mean without removing outliers: 100.40179124497
## Mean if we remove outliers: 100.369856402439
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 1 from 32800 observations
## Proportion (%) of outliers: 0.00304878048780488
## Mean of the outliers: 361.84
## Mean without removing outliers: 100.369856402439
## Mean if we remove outliers: 100.361884508674
## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 0 from 32799 observations
## Proportion (%) of outliers: 0
## Mean of the outliers: NaN
## Mean without removing outliers: 100.361884508674
## Mean if we remove outliers: 100.361884508674
## Outliers successfully removed



plot(DHI$`APORTACION (m3/s)`,
     xlab = "Date",
     ylab= "Aportacion [m³/s]",
     type = "l")



Aplicando esta funcion a los datos de aportación vemos como se han eliminado valores… pero observando los datos… nos damos cuenta de que nuestros datos siguen teniendo demasiado “ruido”. Antes de proseguir hay que rellenar los huecos que hemos ido generando.



2.1.2 Rellenando huecos



Tras hacer la limpieza de los datos, se han generado gran cantidad de NA’s dentro de nuestro dataset.



sum(is.na(DHI$`APORTACION (m3/s)`))
## [1] 7180



A continuación añadimos datos interpolando. Con el paquete imputeTS se puede añadir los datos que faltan interpolando de maneras diferentes:

  • linear (por defecto)
  • spline
  • stine



library(imputeTS)

#Rellenamos Huecos
DHI$`APORTACION (m3/s)`<-  na.interpolation(DHI$`APORTACION (m3/s)`)


#Comprobamos
sum(is.na(DHI$`APORTACION (m3/s)`))
## [1] 0



2.1.3 SMA sobre los datos 15-minutales



Para eliminarlo aplicamos una moving average. Tenemos en cuenta de que tenemos un total de 96 datos diarios. probamos la Moving Average para varios valores.



library(TTR)

for(n in c(1,seq(48,192,length.out = 4))){
  Mavg<- SMA(DHI$`APORTACION (m3/s)`, n) 
  
  plot(DHI$DATE, Mavg,
       xlab = "Date",
       ylab= "Aportacion [m³/s]", 
       type = "l",
       main = paste0("Aportación en Belesar. MA= ", n))
}



Se puede apreciar como de este mode se elimina bastante ruido de la señal de aportación. Pero no sé hasta que punto estamos respetando los datos… Para continuar se decide que utilizaremos una Moving Average con 96 puntos(cantidad de datos generados en 1 dia)… que sería algo parecido a una media diaria (aunque no es lo mismo).



2.1.4 Limpiando datos de nivel del embalse



Tras ejecutar la funcion en buscar de outliers vemos que no hay que en este dataset hay menos datos erroneos.



outlierKD(DHI,`NIVEL EMBALSE (msnm)`)
## Outliers identified: 1 from 39978 observations
## Proportion (%) of outliers: 0.00250137575666617
## Mean of the outliers: 0
## Mean without removing outliers: 310.775323928161
## Mean if we remove outliers: 310.783097781224
## Outliers successfully removed

DHI$`NIVEL EMBALSE (msnm)`<-  na.interpolation(DHI$`NIVEL EMBALSE (msnm)`)



2.1.5 Correlación



A continuación comprobaremos la correlacion que existe entre la aportacion y la variacion del nivel de la presa de Belesar.



x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(c(0,diff(DHI$`NIVEL EMBALSE (msnm)`)*8000+150), x)

aportacion_SMA<- aportacion_SMA[!is.na(aportacion_SMA)]
nivel_SMA<- nivel_SMA[!is.na(nivel_SMA)]


plot(aportacion_SMA,type = "l",
     xlab =" " ,ylab = "",xaxt= 'n',yaxt='n',
     main=paste0("Variacion del nivel y aportacion\n Correlación de: ", 
                 round(cor(aportacion_SMA, nivel_SMA), 3)))
lines(nivel_SMA, col = "red")



2.1.6 Cross-Correlation



Otro análisis previo que se puede hacer es el de la crosscorrelation… para observar con que desfase se obtiene la mejor correlation. Esto es super útil para analizar correlacione de variables que estan desfasadas en el tiempo.



ccf_belesar<- ccf(aportacion_SMA, nivel_SMA, lag.max = 5000)



#Máxima correlación
max(ccf_belesar$acf)
## [1] 0.5892262
#Con cuanto desfase se produce la máxima correlación. 
ccf_belesar$lag[which.max(ccf_belesar$acf)]
## [1] 34



De esta manera podemos darnos cuenta de que estas variables obtienen la mejor correlacion (0.58) para un retraso en la aportación de 34 valores, teniendo en cuenta que nuestro dataset es de datos 15 minutales esto supone 9 horas, es decir, los datos de aportación influyen en el nivel de la presa con un retardo de medio día.



2.1.7 Obtención de medias horarias



A continuación construiremos una tabla con las variables aportacion y nivel obtenidas por diferentes métodos:



  • Haciendo SMA sobre el dataset completo (Hecho anteriormente)
  • Haciendo medias horarias
  • Haciendo SMA sobre las medias horarias



#Retrasamos hacia detrás todos los datos de lluvia porque la acumulada de toda la hora la suma en la hora siguiente....
DHI$`LLUVIA ACUMULADA DÍA (l/m2)`<- lead(DHI$`LLUVIA ACUMULADA DÍA (l/m2)`)

Desacumular_lluvia_2018<- DHI[which(year(DHI$DATE)==2018),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
                                                                         lluvia=ifelse(desacumulada>=0, desacumulada, 0))


Desacumular_lluvia_2019<- DHI[which(year(DHI$DATE)==2019),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
                                                                              lluvia=ifelse(desacumulada>=0, desacumulada, 0))
Desacumular_lluvia<- rbind(Desacumular_lluvia_2018, Desacumular_lluvia_2019)


Medias_horarias_2018<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2018,]) %>% group_by(yday(DATE), hour(DATE))  %>% 
  summarize(., Acum_horaria=sum(lluvia, na.rm = T),
            aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
            nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))

Medias_horarias_2019<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2019,]) %>% group_by(yday(DATE),hour(DATE))  %>% 
  summarize(., Acum_horaria=sum(lluvia, na.rm = T),
            aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
            nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))

Medias_horarias<- rbind(Medias_horarias_2018, Medias_horarias_2019)
vector_Date<- seq(range(DHI$DATE)[1],range(DHI$DATE)[2],
                  by="hour")

Lluvia_acum_horaria<- as.data.frame(cbind(as.character(vector_Date[2:length(vector_Date)]), 
                                          Medias_horarias[,3:5]))



Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`<- ymd_hms(Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`)


colnames(Lluvia_acum_horaria)<- c("Date", "Lluvia_mm", "aport_mean", "nivel_mean")



2.1.7.1 Comprobación mejor SMA sobre la media horaria



Comprobacoin de máxima correlacion desplazada para diferentes valores de SMA.



for (n in seq(6,24*3, by=6)) {
  y<- n
  aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
  nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)
  
  
  
  Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria, 
                                  
                                  aportacion_mean_SMA,
                                  nivel_mean_SMA))
  
  
  
  colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_mean_SMA",  "nivel_mean_SMA")
  Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]
  
  
  
  ccf_aport_mean<- ccf(Tabla_DHI$aport_mean_SMA,
                       Tabla_DHI$nivel_mean_SMA, 
                       lag.max = 5000, 
                       plot = T)
  
  print(paste0("Máxima correlacion de ",
               round(max(ccf_aport_mean$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_aport_mean$lag[which.max(ccf_aport_mean$acf)]/24, digits = 2), 
               " días. Con un SMA de :", 
               n, 
               " horas"))
  
}

## [1] "Máxima correlacion de 0.609 .Para un desfase de: -76.79 días. Con un SMA de :6 horas"

## [1] "Máxima correlacion de 0.616 .Para un desfase de: -76.58 días. Con un SMA de :12 horas"

## [1] "Máxima correlacion de 0.621 .Para un desfase de: -76.42 días. Con un SMA de :18 horas"

## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días. Con un SMA de :24 horas"

## [1] "Máxima correlacion de 0.627 .Para un desfase de: -76.29 días. Con un SMA de :30 horas"

## [1] "Máxima correlacion de 0.629 .Para un desfase de: -76.25 días. Con un SMA de :36 horas"

## [1] "Máxima correlacion de 0.632 .Para un desfase de: -76.21 días. Con un SMA de :42 horas"

## [1] "Máxima correlacion de 0.634 .Para un desfase de: -76.08 días. Con un SMA de :48 horas"

## [1] "Máxima correlacion de 0.636 .Para un desfase de: -75.92 días. Con un SMA de :54 horas"

## [1] "Máxima correlacion de 0.638 .Para un desfase de: -75.88 días. Con un SMA de :60 horas"

## [1] "Máxima correlacion de 0.64 .Para un desfase de: -75.79 días. Con un SMA de :66 horas"

## [1] "Máxima correlacion de 0.642 .Para un desfase de: -75.67 días. Con un SMA de :72 horas"



Se puede observar como casi todos los casos probando diferentes SMA’s llegan a la conclusión de que la correlación máxima se tiene para 3 días y poco entre aportación y nivel… Se logra mejorar la correlacion aumentando el SMA. La mejor correlación aparece para un desfase de entorno a 75 días. Lo cual no tiene mucho sentido. o por lo menos a nosotros no nos vale…



2.1.8 Creamos Dataset final



Añadimos ahora los SMA’s de todo el dataset y el SMA sobre la media horaria.



#se selecciona 96 porque es la cantidad de datos que corresponden a 1 día. 
x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(DHI$`NIVEL EMBALSE (msnm)`, x)

aportacion_horaria<- aportacion_SMA[DHI$DATE%in%vector_Date]
nivel_horario<- nivel_SMA[DHI$DATE%in%vector_Date]

#Se ponen 24 horas para que equivalgan a 1 día 
y<- 24
aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)





Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria, 
                                aportacion_horaria[2:length(aportacion_horaria)],
                                nivel_horario[2:length(nivel_horario)], 
                                aportacion_mean_SMA,
                                nivel_mean_SMA))



colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_SMA", "nivel_SMA", 
                        "aport_mean_SMA",  "nivel_mean_SMA")

Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]



2.1.9 Ploteos de aportacion y nivel.



A continuación representamos la aportacion y el nivel. para los tres casos



2.1.9.1 SMA sobre media horaria



 k<- max(Tabla_DHI$nivel_mean_SMA)/max(Tabla_DHI$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_mean_SMA/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_mean_SMA),
                                                        max(Tabla_DHI$nivel_mean_SMA),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean_SMA),
                                    max(Tabla_DHI$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Nivel y aportación. SMA (24 horas) de la media horaria")

ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA, 
    Tabla_DHI$nivel_mean_SMA, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.583 .Para un desfase de: -41.67 días"



2.1.9.2 SMA sobre datos 15-minutales



 k<- max(Tabla_DHI$nivel_SMA)/max(Tabla_DHI$aport_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_SMA/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_SMA),
                                                        max(Tabla_DHI$nivel_SMA),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_SMA),
                                    max(Tabla_DHI$aport_SMA),
                                    by=20)) + 
    ggtitle("Nivel y aportación. SMA (96 datos = 1 día) de los datos 15-minutales")

ccf_DHI<- ccf(Tabla_DHI$aport_SMA, 
    Tabla_DHI$nivel_SMA, 
    lag.max = 4000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días"



2.1.9.3 Media horaria sin SMA



 k<- max(Tabla_DHI$nivel_mean)/max(Tabla_DHI$aport_mean)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_mean/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_mean),
                                                        max(Tabla_DHI$nivel_mean),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean),
                                    max(Tabla_DHI$aport_mean),
                                    by=20)) + 
    ggtitle("Nivel y aportación. Medias horarias")

ccf_DHI<- ccf(Tabla_DHI$aport_mean, 
    Tabla_DHI$nivel_mean, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.54 .Para un desfase de: -41.67 días"



Por lo anteriormente expuesto… se llega a la conclusión de que el mejor método para relacionar la aportación y el nivel es haciendo una moving average sobre las medias horarias… un SMA de 142 valores… entorno a 6 dias.



2.1.10 Correlación entre la lluvia y la aportación observada



Otra comprobacion que considero pertinente es ver si la lluvia y la aportacion correlan debidamente…



 k<- max(Tabla_DHI$Lluvia_mm)/max(Tabla_DHI$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Lluvia Horaria [l/m²]", 
                                           breaks = seq(min(Tabla_DHI$Lluvia_mm),
                                                        max(Tabla_DHI$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean_SMA),
                                    max(Tabla_DHI$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion")

ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA, 
    Tabla_DHI$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.109 .Para un desfase de: 1.71 días"



Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12), ]
 k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI_cut, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
                                                        max(Tabla_DHI_cut$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
                                    max(Tabla_DHI_cut$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion. Diciembre")

ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA, 
    Tabla_DHI_cut$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.371 .Para un desfase de: 1.62 días"



Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12 | year(Tabla_DHI$Date)==2019), ]
 k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI_cut, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
                                                        max(Tabla_DHI_cut$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
                                    max(Tabla_DHI_cut$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion. Diciembre-Enero-Febrero")

ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA, 
    Tabla_DHI_cut$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))
## [1] "Máxima correlacion de 0.332 .Para un desfase de: 1.71 días"

Se puede observar como la correlación entre la aportación y la lluvia es apriori muy mala… habrá que hacer trikiñuelas para poder salvar esto… En principio me parece normal que estas variables tengan mala correlación… la aportación del río si que dependerá e la lluvia… pero influyen muchos más factores… además. La lluvia registrada en belesar no tiene porqué ser la mejor para estimar la aportación del río que llega a belesar… puede estar lluviendo muchos kilometros hacia arriba en el valle…



Antes de proseguir guardamos este histórico que hemos construido.



path_dataDHI <- here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS')

#saveRDS(Tabla_DHI,path_dataDHI)



2.2 Datos WRF



Los datos del histórico de WRF están organizados en formato lista… hay para cada localización 2 dataframes;

  • D1: WRF predicción del mismo día
  • D2: WRF predicción de mañana



A continuación se presenta el loop y las funciones que se emplean para generar todo el histórico WRF de Belesar… Cabe mencionar que contamos con una carpeta con Archivos RDS ordenados por fechas para la localización de belesar.



Usamos las siguientes funciones.



lon_lat_df_ls<- function(parque_list){
  fecha<- names(parque_list)
  fecha[!str_detect(fecha, " ")]<- paste0(fecha[!str_detect(fecha, " ")], " 00:00:00")
  
  n_lon<- parque_list[[1]]$lon
  n_lat<- parque_list[[1]]$lat
  
  nombres<- paste0(n_lon,"__",n_lat)
  
  
  list_localizaciones<- list()
  for (localizaciones in 1:length(parque_list[[1]][,1]) ) {
    l<- lapply(parque_list, function(x) x[localizaciones,])
    df <- data.frame(matrix(unlist(l), nrow=length(l), byrow=T))
    df<- as.data.frame(cbind(ymd_hms(fecha), df))
    colnames(df)<-c("fechas", colnames(parque_list[[1]]))
    list_localizaciones[[localizaciones]]<- df
  }
  
  names(list_localizaciones)<- nombres
  return(list_localizaciones)
  
  
}
uv_transformation<- function(tabla_comp){
  
  nombres<- colnames(tabla_comp)
  
  u10<- tabla_comp$U10_MEAN
  v10<- tabla_comp$V10_MEAN
  
  u10_max<- tabla_comp$U10_MAX
  v10_max<- tabla_comp$V10_MAX
  
  wind_abs <- sqrt(u10^2 + v10^2)
  wind_dir_rad <-  atan2(u10/wind_abs, v10/wind_abs) 
  wind_dir_deg1 <-  wind_dir_rad * 180/pi 
  wind_dir_deg2 <-  wind_dir_deg1+ 180 
  
  
  wind_abs_max <- sqrt(u10_max^2 + v10_max^2)
  wind_dir_rad_max <-  atan2(u10_max/wind_abs_max, v10_max/wind_abs_max) 
  wind_dir_deg1_max <-  wind_dir_rad_max * 180/pi 
  wind_dir_deg2_max <-  wind_dir_deg1_max+ 180 
  
  tabla_comp<- as.data.frame(cbind(tabla_comp,wind_abs,wind_dir_deg2,
                                   wind_abs_max,wind_dir_deg2_max))
  colnames(tabla_comp)<- c(nombres,"WS","WD","WS_MAX","WD_MAX")
  tabla_comp$U10_MAX<- NULL
  tabla_comp$V10_MAX<- NULL
  tabla_comp$U10_MEAN<- NULL
  tabla_comp$V10_MEAN<- NULL
  return(tabla_comp)
  
}
extract_rain_data<- function(Belesar_lolat_df){
  
  rain<- Belesar_lolat_df[,c(1,2,3,4,5)]
  rain$pre_acum<- rain$RAINC+rain$RAINNC
  rain$RAINC<- NULL
  rain$RAINNC<- NULL
  
  prep_hourly<- vector()
  for (i in 1:length(rain$pre_acum)) {
    if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
      prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
    }
    
  }
  rain$prep_hourly<- prep_hourly 
  
  return(rain)
}
extract_rain_data2<- function(Belesar_lolat_df){
  
  rain<- Belesar_lolat_df[,c(1,2,3,4,5,6,10,17,20)]
  rain$pre_acum<- rain$RAINC+rain$RAINNC
  
  prep_hourly<- vector()
  for (i in 1:length(rain$pre_acum)) {
    if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
      prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
    }
    
  }
  rain$prep_hourly<- prep_hourly 
  
  return(rain)
}
Return_periodo_Belesar<- function(){
  All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'), 
                                 recursive = F, full.names = T)
  RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]
  
  RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]
  
  
  Belesar_data<- readRDS(RDS_Belesar1[1])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
  
  fecha_ini<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[1]
  
  Belesar_data<- readRDS(RDS_Belesar1[length(RDS_Belesar)])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
  
  fecha_last<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[length(Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas)]
  
  periodo_WRF<- seq(fecha_ini, fecha_last, by="hour")
  
  Tabla_WRF<- as.data.frame(matrix(ncol = 10, nrow = length(periodo_WRF)))
  colnames(Tabla_WRF)<- colnames(Belesar_rain[[1]])
  Tabla_WRF$fechas<- periodo_WRF
  return(Tabla_WRF)
}



#Listamos archivos dentro de la carpeta de Belesar
All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'),
                               recursive = F, full.names = T)


#Detectamos cuales son RDS
RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]

#Eliminamos los RDS que no son de WRF
RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]

#Construimos Lista para cada instante de tiempo
Lista_localizacion<- list() 
for (i in 1:length(RDS_Belesar1)) {
  Belesar_data<- readRDS(RDS_Belesar1[i])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data2)
  Lista_localizacion[[i]]<- Belesar_rain
}

#Nombramos la lista
names_fechas<- sapply(RDS_Belesar1, function(x){
  r<- str_split(x, "/")
  return(str_remove(str_remove(r[[1]][length(r[[1]])], ".RDS"), "Belesar_"))
})
names(Lista_localizacion)<- names_fechas


#Creamos una lista por localización
Lista_localizacion2<- list()
for (i in 1:length(Lista_localizacion[[1]])) {
  Lista_localizacion2[[i]]<- lapply(Lista_localizacion, 
                               function(x) return(x[[i]]))
}

#Nombramos la lista
names(Lista_localizacion2)<- names(Lista_localizacion[[1]])

#Guardamos la lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
#saveRDS(Lista_localizacion2, nombre_lista)

#Cargamos lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
Lista_total1<- readRDS(nombre_lista)

#Juntamos todos los Dataframes
Lista_total_MF<- lapply(Lista_total1, function(x) bind_rows(x))

#creamos dos data.frames... uno para D1 y otro para D2
Lista_d1_d2_loc<- list()
for (i in 1:length(Lista_total_MF)) {
  p<- Lista_total_MF[[i]]
  d1<- p[duplicated(p$fechas),]
  d1<-d1[!duplicated(d1$fechas),]
  d2<- p[!duplicated(p$fechas),]
  
  d2_qneed1<-d2[!(d2$fechas%in%d1$fechas),]
  
  
  d1_2<-bind_rows(d1,d2_qneed1)
  d1_2<-d1_2[order(d1_2$fechas),]
  
  d2<-d2[order(d2$fechas),]
  
  d1_2$pre_acum<- NULL
  d2$pre_acum<- NULL
  
  colnames(d1_2)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
  
  colnames(d2)<-  c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
  
  lista_loc_d12<- list(d1_2,d2)
  names(lista_loc_d12)<- c("D1", "D2")
  
  Lista_d1_d2_loc[[i]]<- lista_loc_d12
}

#nombramos la lsta
names(Lista_d1_d2_loc)<- names(Lista_total_MF)

Tabla_periodo1<- Return_periodo_Belesar()
colnames(Tabla_periodo1)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")

Lista_d1_d2_loc2<- list()
for (j in 1:length(Lista_d1_d2_loc)) {
  prueba_list<- Lista_d1_d2_loc[[j]]
  lista_retorno<- list()
  for(i in 1:2){
    prueba<- prueba_list[[i]]
    Tabla_periodo<- Tabla_periodo1
    Tabla_periodo$LON[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LON
    Tabla_periodo$LAT[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LAT
    Tabla_periodo$RAINC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINC
    Tabla_periodo$RAINNC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINNC
    Tabla_periodo$RAINSH[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINSH
    Tabla_periodo$T02_MEAN[match(prueba$Date,Tabla_periodo$Date)] <- prueba$T02_MEAN
    Tabla_periodo$PSFC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$PSFC
    Tabla_periodo$WS_MAX[match(prueba$Date,Tabla_periodo$Date)] <- prueba$WS_MAX
    Tabla_periodo$prep_hourly[match(prueba$Date,Tabla_periodo$Date)] <- prueba$prep_hourly
    lista_retorno[[i]]<- Tabla_periodo
  }
  names(lista_retorno)<- c("D1", "D2")
  Lista_d1_d2_loc2[[j]]<- lista_retorno 
  
  
}
names(Lista_d1_d2_loc2)<- names(Lista_d1_d2_loc)

#Creamos lista con las variables afinadas
Lista_d1_d2_loc3<- lapply(Lista_d1_d2_loc2, function(x){
  x[[1]]$RAINC<- NULL
  x[[1]]$RAINNC<- NULL
  x[[1]]$RAINSH<- NULL
  x[[2]]$RAINC<- NULL
  x[[2]]$RAINNC<- NULL
  x[[2]]$RAINSH<- NULL
  
  return(x)})
  
names(Lista_d1_d2_loc3)<- names(Lista_d1_d2_loc2)

#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS')
#saveRDS(Lista_d1_d2_loc3,path_hist_WRF)



2.3 Juntando históricos



Juntamos a cada data.frame de WRF la tabla afinada de DHI.



#Usamos dplyr para juntar ambos datasets

Belesar_WRF<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS'))
Belesar_DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS'))

df2<- Belesar_DHI

Belesar_Merge<- list()
for (j in 1:length(Belesar_WRF)) {
  lista_retorno<- list()
  for(i in 1:2){
    df1<-  Belesar_WRF[[j]][[i]]
    Merge_table<- left_join(df1, df2, by=c("Date"))
    lista_retorno[[i]]<- Merge_table
  }
  names(lista_retorno)<- c("D1", "D2")
  Belesar_Merge[[j]]<- lista_retorno 
}
names(Belesar_Merge)<- names(Belesar_WRF)
#Belesar merge completecases
Belesar_Merge_cc<- list()
for (j in 1:length(Belesar_Merge)) {
    lista_retorno<- list()
  for(i in 1:2){
    df1<-  Belesar_Merge[[j]][[i]]
    df1$Acumulated<- NULL
    Table_fine<- df1[complete.cases(df1),]
    lista_retorno[[i]]<- Table_fine
  }
  names(lista_retorno)<- c("D1", "D2")
  Belesar_Merge_cc[[j]]<- lista_retorno 
}
names(Belesar_Merge_cc)<- names(Belesar_Merge)



#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
#saveRDS(Belesar_Merge_cc,path_hist_WRF)



3 Comprobacion de correlacion aplicando diferentes regresiones lineales



#importar
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)


#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data, function(x){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/10/01")
    fecha_end<- ymd("2019/02/01")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
})
cut_predict<- lapply(clean_data, function(x){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2019/02/01")
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
})



lista_TL<- list()
for (i in 1:length(cut_train)) {
  lista_d1d2<- list()
  for (j in 1:2) {
    data_predict<- cut_train[[i]][[j]]
    data_predict2<- cut_predict[[i]][[j]]
    
    fit_1  <- lm(Lluvia_mm  ~ prep_hourly, data = data_predict)
    fit_2  <- lm(Lluvia_mm  ~ prep_hourly + T02_MEAN, data = data_predict)
    fit_3  <- lm(Lluvia_mm  ~ prep_hourly + PSFC, data = data_predict)
    fit_4  <- lm(Lluvia_mm  ~ prep_hourly + WS_MAX, data = data_predict)
    fit_5  <- lm(Lluvia_mm  ~ prep_hourly * T02_MEAN, data = data_predict)
    fit_6  <- lm(Lluvia_mm  ~ prep_hourly * PSFC, data = data_predict)
    fit_7  <- lm(Lluvia_mm  ~ prep_hourly * WS_MAX, data = data_predict)
    
    uncorrected<-data_predict2$prep_hourly
    prediction_rain<- predict(fit_1, data.frame(prep_hourly =data_predict2$prep_hourly))
    prediction_rain2<- predict(fit_2, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 T02_MEAN =data_predict2$T02_MEAN))
    prediction_rain3<- predict(fit_3, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 PSFC =data_predict2$PSFC))
    prediction_rain4<- predict(fit_4, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 WS_MAX =data_predict2$WS_MAX))
    prediction_rain5<- predict(fit_5, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 T02_MEAN =data_predict2$T02_MEAN))
    prediction_rain6<- predict(fit_6, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 PSFC =data_predict2$PSFC))
    prediction_rain7<- predict(fit_7, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 WS_MAX =data_predict2$WS_MAX))
    
    observed_rain<-data_predict2$Lluvia_mm 
    
    
    
    tabla_cor<-cbind(cor(uncorrected, observed_rain),
                     cor(prediction_rain, observed_rain),
                     cor(prediction_rain2, observed_rain),
                     cor(prediction_rain3, observed_rain),
                     cor(prediction_rain4, observed_rain),
                     cor(prediction_rain5, observed_rain),
                     cor(prediction_rain6, observed_rain),
                     cor(prediction_rain7, observed_rain))
    
    fit_1  <- svm(Lluvia_mm  ~ prep_hourly, data = data_predict)
    fit_2  <- svm(Lluvia_mm  ~ prep_hourly + T02_MEAN, data = data_predict)
    fit_3  <- svm(Lluvia_mm  ~ prep_hourly + PSFC, data = data_predict)
    fit_4  <- svm(Lluvia_mm  ~ prep_hourly + WS_MAX, data = data_predict)
    fit_5  <- svm(Lluvia_mm  ~ prep_hourly * T02_MEAN, data = data_predict)
    fit_6  <- svm(Lluvia_mm  ~ prep_hourly * PSFC, data = data_predict)
    fit_7  <- svm(Lluvia_mm  ~ prep_hourly * WS_MAX, data = data_predict)
    
    uncorrected<-data_predict2$prep_hourly
    prediction_rain<- predict(fit_1, data.frame(prep_hourly =data_predict2$prep_hourly))
    prediction_rain2<- predict(fit_2, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 T02_MEAN =data_predict2$T02_MEAN))
    prediction_rain3<- predict(fit_3, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 PSFC =data_predict2$PSFC))
    prediction_rain4<- predict(fit_4, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 WS_MAX =data_predict2$WS_MAX))
    prediction_rain5<- predict(fit_5, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 T02_MEAN =data_predict2$T02_MEAN))
    prediction_rain6<- predict(fit_6, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 PSFC =data_predict2$PSFC))
    prediction_rain7<- predict(fit_7, data.frame(prep_hourly =data_predict2$prep_hourly,
                                                 WS_MAX =data_predict2$WS_MAX))
    
    observed_rain<-data_predict2$Lluvia_mm 
    
    tabla_cor<-cbind(tabla_cor,
                     cor(prediction_rain, observed_rain),
                     cor(prediction_rain2, observed_rain),
                     cor(prediction_rain3, observed_rain),
                     cor(prediction_rain4, observed_rain),
                     cor(prediction_rain5, observed_rain),
                     cor(prediction_rain6, observed_rain),
                     cor(prediction_rain7, observed_rain))
    
    colnames(tabla_cor)<- c("uncorrected", "1","2", "3","4","5","6","7",
                            "svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
    lista_d1d2[[j]]<- tabla_cor
  }
  
  names(lista_d1d2)<- c("D1", "D2")
  lista_TL[[i]]<- lista_d1d2
}


colnames(tabla_cor)<- c("uncorrected", "1","2", "3","4","5","6","7",
                        "svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
names(lista_d1d2)<- c("D1", "D2")


Lista_TL_names<-lapply(lista_TL, function(x){
  y<- lapply(x, function(r) {
    r<- as.data.frame(r)
    names(r)<- c("uncorrected", "1","2", "3","4","5","6","7",
                 "svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
    return(r)
  })
  names(y)<- c("D1", "D2")
  return(y)
})
names(Lista_TL_names)<- names(cut_predict)


x<- lapply(Lista_TL_names, function(x){bind_rows(as.data.frame(x))})
y<- bind_rows(x, .id="id")

cor_table<- as.data.frame(matrix(ncol = 3, nrow = 43))
colnames(cor_table)<- c("id","corr", "name")

for (i in 2:length(y[1,])) {
  f<- as.data.frame(cbind(y[which.max(y[,i]), c(1,i)], names( y[which.max(y[,i]), c(1,i)])[2]))
  colnames(f)<- c("id","corr", "name")
  cor_table<- rbind(cor_table,f )
}

cor_table<- cor_table[complete.cases(cor_table),]

punto_Belesar<- c(42.628577,-7.713948)
extract_lat_lon<- lapply(str_split(cor_table$id,"__"), function(x){
  return(as.data.frame(cbind(x[[1]],x[[2]])))
})

extract_lat_lon<- as.data.frame(bind_rows(extract_lat_lon))


cor_table<- as.data.frame(cbind(extract_lat_lon, cor_table[,2:3]))

colnames(cor_table)<- c("LON","LAT", "Corr", "Method")
#Librari Geosphere para medir la distancia entre dos puntos a partir de sus coordenadas
library(geosphere)
vec_dist<- vector()
for (i in 1:length(cor_table[,1])) {
  vec_dist[i]<- distm(c(-7.713948, 42.628577), c(as.numeric(cor_table[i,1]), 
                                                 as.numeric(cor_table[i,2])), 
                      fun = distHaversine)
  
}

vec_dist<- round(vec_dist/1000, digits = 1)

cor_table<- as.data.frame(cbind(cor_table, vec_dist))
cor_table$LON<-round(as.numeric(cor_table$LON), digits = 2)
cor_table$LAT<-round(as.numeric(cor_table$LAT), digits = 2)
cor_table$Corr<-round(as.numeric(cor_table$Corr), digits = 3)

cor_table



De esta manera llegamos a la conclusión de que empleando diferentes variables entregadas por WRF, tales como presión superficial, viento, y temperatura. No existe una regresión lineal o SVM que pueda mejorar la predicción de la precipitación ofrecida por WRF.

4 Correlación lluvia WRF y lluvia observada

Pasamos ahora a analizar la lluvia ofrecida por WRF comparada con la registrada en Belesar. Ya hemos visto que la correlación mayor se obtiene para un punto a 23 km de distancia.

4.1 Octubre-Febrero

#Empleamos los datos a partir de octubre y eliminamos datos negativos de 
# precipitacion WRF que hemos visto que hay alguno. 
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) {
    r<- r[which(r$Date > ymd("2018/10/01")),]
    r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly) 
    return(r)
  })
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_oct)) {
  Corr<- cor(clean_data_oct[[i]][[1]]$prep_hourly, 
      clean_data_oct[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
                      c(as.numeric(LON),
                        as.numeric(LAT)),
                      fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order
#Plot_rain3 necesita columna Date y columna rain
plot_rain3<- function(data,data2, titulo){
  ggplot(data=data, aes(x=Date))+
    geom_line(aes(y=rain), stat="identity")+
    xlab("Date")+ylab("Lluvia por hora [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank()) +
    geom_line(data= data2, aes(x=Date, y=rain), color="red", alpha=0.5)+
    geom_text(aes(as.POSIXct(ymd("2019/02/01")),10, label=paste0("Cor: ",round(cor(data$rain,data2$rain), 
                                                                               digits = 2)))) +
    ggtitle(titulo) + 
     theme(plot.title = element_text(hjust = 0.5))
  
}
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
                              clean_data_oct[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
                              clean_data_oct[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_oct)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

4.2 Enero-Febrero

Es reseñable que disminuyendo el tiempo de ploteo la correlación aumenta.

clean_data_jan<-lapply(clean_data_oct, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/01")),])
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
  Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly, 
             clean_data_jan[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}

Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]

Cor_order
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

4.3 15 de Enero- 10 de febrero

Apuramos más el rango para que únicamente coja un periodo lluvioso como fue el de el 15 de enero al pocos de febrero

clean_data_jan<-lapply(clean_data_oct, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/15") & r$Date < ymd("2019/02/10")),])
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
  Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly, 
             clean_data_jan[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

5 Comprobación de Medias Diarias

A continuacion se realizan las medias diarias de todos los datos y se comprueba la correlación entre los datos de lluvia y aportación.

Daiyli_avgs<- lapply(clean_data_oct, function(x){
  x[[1]] %>% group_by(year(Date), yday(Date)) %>%  
    summarize(rWRF_acum= sum(prep_hourly),
              rDHI_acum= sum(Lluvia_mm),
              rWRF_mean= sum(prep_hourly),
              rDHI_mean= sum(Lluvia_mm),
              aport_SMA_diario=mean(aport_SMA),
              aport_mean_diario=mean(aport_mean),
              aport_mean_SMA_diario=mean(aport_mean_SMA),
              nivel_SMA_diario=mean(nivel_SMA),
              nivel_mean_diario=mean(nivel_mean),
              nivel_mean_SMA_diario=mean(nivel_mean_SMA))

  
  })

Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
  
  x$Date<- ymd(ifelse(x$`year(Date)`==2018,
                      as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
                      as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
  x$`year(Date)`<- NULL
  x$`yday(Date)`<- NULL
  return(x)
})

Cor_rain_place<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  Corr1<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
             Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr2<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr3<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr4<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum)
  
  Corr5<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr6<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr7<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr8<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, round(Corr1, digits = 2),
                                           round(Corr2, digits = 2),
                                           round(Corr3, digits = 2),
                                           round(Corr4, digits = 2),
                                           round(Corr5, digits = 2),
                                           round(Corr6, digits = 2),
                                           round(Corr7, digits = 2),
                                           round(Corr, digits = 2),
                                           Dist))
  
  
}

colnames(Cor_rain_place)<- c("LON","LAT","AWRF_aport1",
                             "AWRF_aport2", "AWRF_aport3",
                             "AWRF_ADHI","MWRF_aport1",
                             "MWRF_aport2", "MWRF_aport3",
                             "MWRF_MDHI","Dist")

Se aprecia que los datos de aportación son muy malos

Cor_order<- Cor_rain_place[order(Cor_rain_place$Dist),]
Cor_order
Cor_rain_place_ccf<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  ccfr1<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr2<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr3<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr4<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum, plot = FALSE)
  
  ccfr5<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr6<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr7<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr8<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean, plot = FALSE)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place_ccf[i,]<- as.data.frame(cbind(LON, LAT, 
                                               round(max(ccfr1$acf), digits = 2),
                                               round(max(ccfr2$acf), digits = 2),
                                               round(max(ccfr3$acf), digits = 2),
                                               round(max(ccfr4$acf), digits = 2),
                                               round(max(ccfr5$acf), digits = 2),
                                               round(max(ccfr6$acf), digits = 2),
                                               round(max(ccfr7$acf), digits = 2),
                                               round(max(ccfr8$acf), digits = 2), Dist))
  print(paste("Mejor corr para un retardo de: ", 
              ccfr1$lag[which.max(ccfr1$acf)],
              ccfr2$lag[which.max(ccfr2$acf)],
              ccfr3$lag[which.max(ccfr3$acf)],
              ccfr4$lag[which.max(ccfr4$acf)],
              ccfr5$lag[which.max(ccfr5$acf)],
              ccfr6$lag[which.max(ccfr6$acf)],
              ccfr7$lag[which.max(ccfr7$acf)],
              ccfr8$lag[which.max(ccfr8$acf)], sep = " "))  
  
}
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -6 -5 -6 0 -6 -5 -6 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -4 -3 0 -3 -4 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -13 -12 -13 0 -13 -12 -13 0"
colnames(Cor_rain_place_ccf)<-  c("LON","LAT","AWRF_aport1",
                                  "AWRF_aport2", "AWRF_aport3",
                                  "AWRF_ADHI","MWRF_aport1",
                                  "MWRF_aport2","MWRF_aport3",
                                  "MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place_ccf[order(Cor_rain_place_ccf$Dist),]
Cor_order

5.1 Reducimos el tiempo

Antes estabamos usando los datos desde octubre hasta febrero… Ahora empleamos solo los datos desde diciembre hasta febrero.

clean_data_dec<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2018/12/01")),])
})
Daiyli_avgs<- lapply(clean_data_dec, function(x){
  x[[1]] %>% group_by(year(Date), yday(Date)) %>%  
    summarize(rWRF_acum= sum(prep_hourly),
              rDHI_acum= sum(Lluvia_mm),
              rWRF_mean= sum(prep_hourly),
              rDHI_mean= sum(Lluvia_mm),
              aport_SMA_diario=mean(aport_SMA),
              aport_mean_diario=mean(aport_mean),
              aport_mean_SMA_diario=mean(aport_mean_SMA),
              nivel_SMA_diario=mean(nivel_SMA),
              nivel_mean_diario=mean(nivel_mean),
              nivel_mean_SMA_diario=mean(nivel_mean_SMA))

  
  })

Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
  
  x$Date<- ymd(ifelse(x$`year(Date)`==2018,
                      as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
                      as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
  x$`year(Date)`<- NULL
  x$`yday(Date)`<- NULL
  return(x)
})

Cor_rain_place<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  Corr1<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
             Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr2<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr3<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr4<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum)
  
  Corr5<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr6<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr7<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr8<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, round(Corr1, digits = 2),
                                           round(Corr2, digits = 2),
                                           round(Corr3, digits = 2),
                                           round(Corr4, digits = 2),
                                           round(Corr5, digits = 2),
                                           round(Corr6, digits = 2),
                                           round(Corr7, digits = 2),
                                           round(Corr, digits = 2),
                                           Dist))
  
  
}

colnames(Cor_rain_place)<- c("LON","LAT","AWRF_aport1",
                             "AWRF_aport2", "AWRF_aport3",
                             "AWRF_ADHI","MWRF_aport1",
                             "MWRF_aport2", "MWRF_aport3",
                             "MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place[order(Cor_rain_place$Dist),]
Cor_order
Cor_rain_place_ccf<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  ccfr1<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr2<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr3<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr4<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum, plot = FALSE)
  
  ccfr5<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr6<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr7<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr8<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean, plot = FALSE)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place_ccf[i,]<- as.data.frame(cbind(LON, LAT, 
                                               round(max(ccfr1$acf), digits = 2),
                                               round(max(ccfr2$acf), digits = 2),
                                               round(max(ccfr3$acf), digits = 2),
                                               round(max(ccfr4$acf), digits = 2),
                                               round(max(ccfr5$acf), digits = 2),
                                               round(max(ccfr6$acf), digits = 2),
                                               round(max(ccfr7$acf), digits = 2),
                                               round(max(ccfr8$acf), digits = 2), Dist))
  print(paste("Mejor corr para un retardo de: ", 
              ccfr1$lag[which.max(ccfr1$acf)],
              ccfr2$lag[which.max(ccfr2$acf)],
              ccfr3$lag[which.max(ccfr3$acf)],
              ccfr4$lag[which.max(ccfr4$acf)],
              ccfr5$lag[which.max(ccfr5$acf)],
              ccfr6$lag[which.max(ccfr6$acf)],
              ccfr7$lag[which.max(ccfr7$acf)],
              ccfr8$lag[which.max(ccfr8$acf)], sep = " "))  
  
}
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -6 -5 -6 0 -6 -5 -6 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -4 -3 0 -3 -4 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -13 -12 -13 0 -13 -12 -13 0"
colnames(Cor_rain_place_ccf)<-  c("LON","LAT","AWRF_aport1",
                                  "AWRF_aport2", "AWRF_aport3",
                                  "AWRF_ADHI","MWRF_aport1",
                                  "MWRF_aport2","MWRF_aport3",
                                  "MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place_ccf[order(Cor_rain_place_ccf$Dist),]
Cor_order

A continuación extraemos la media y el máximo de cada columna.

#Media de cada columna
colMeans(Cor_rain_place_ccf[3:(length(Cor_rain_place_ccf)-1)])
## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##   0.2883721   0.2939535   0.2890698   0.6567442   0.2883721   0.2939535 
## MWRF_aport3   MWRF_MDHI 
##   0.2890698   0.6567442
#función ColMAX
colMax <- function(data) sapply(data, max, na.rm = TRUE)

colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]
## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##        0.49        0.50        0.49        0.73        0.49        0.50 
## MWRF_aport3   MWRF_MDHI 
##        0.49        0.73

De aquí podemos extraer las siguientes conclusiones:

  • No hay diferencia en la correlación usando la lluvia acumulada o la lluvia media diaria.
  • Para realizar correlación lluvia WRF y aportación el mejor metodo es el 3, es decir, haciendo SMA sobre los datos medios horarios de aportación. A partir de ahora emplearemos AWRF_aport3 y AWRF_ADHI, como indicadores para la correlacion lluvia WRF-aportación y lluviaWRF- lluvia observada
  • Existe un punto con una correlación entre lluvia registrada y lluvia WRF de 92%. BUENAS NOTICIAS. Esto nos da información de que el modelo está haciendo bien su trabajo.
#función ColMAX
which_colMax <- function(data) sapply(data, which.max)

which_colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]
## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##          18          18          18           1          18          18 
## MWRF_aport3   MWRF_MDHI 
##          18           1

Vemos que en cuanto a correlación entre WRF y aportación el mejor punto es el numero 40.

Cor_rain_place_ccf[40,]

Apreciamos que se encuentra a una distancia de 82 km. Esto es peligroso.

En cuanto a la lluvia…

Cor_rain_place_ccf[4,]

Vemos que se encuentra a 52 km.

Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T),c("AWRF_aport3","AWRF_ADHI", "Dist")]

Por lo visto en la tabla anterior, habrá que buscar un punto que tenga buenas correlaciones para ambas comparaciones y además esté relativamente cerca. En mi opinión estas son las mejores opciones. No sé cual elegir.

Cor_rain_place_ccf[c(19,23,24,14),c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]

Construimos una lista solo con los puntos las variables que nos interesan. Hemos ido reduciendo el rango de tiempo viendo como las correlaciones aumentaban, ahora que hemos escogido el mejor punto vamos a cojer los datos desde finales de octubre.

clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) {
    r<- r[which(r$Date > ymd("2018/10/01")),]
    r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly) 
    return(r)
  })
})
Daiyli_avgs<- lapply(clean_data_oct, function(x){
  x[[1]] %>% group_by(year(Date), yday(Date)) %>%  
    summarize(rWRF_acum= sum(prep_hourly),
              rDHI_acum= sum(Lluvia_mm),
              rWRF_mean= sum(prep_hourly),
              rDHI_mean= sum(Lluvia_mm),
              aport_SMA_diario=mean(aport_SMA),
              aport_mean_diario=mean(aport_mean),
              aport_mean_SMA_diario=mean(aport_mean_SMA),
              nivel_SMA_diario=mean(nivel_SMA),
              nivel_mean_diario=mean(nivel_mean),
              nivel_mean_SMA_diario=mean(nivel_mean_SMA))

  
  })

Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
  
  x$Date<- ymd(ifelse(x$`year(Date)`==2018,
                      as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
                      as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
  x$`year(Date)`<- NULL
  x$`yday(Date)`<- NULL
  return(x)
})

Daiyli_avgs3<-lapply(Daiyli_avgs2[c(19,23,24,14)], function(x)x[, c("Date","rWRF_acum",
                                      "aport_mean_SMA_diario", 
                                      "rDHI_acum")])


head(Daiyli_avgs3[[1]])
range(Daiyli_avgs3[[1]]$Date)
## [1] "2018-10-26" "2019-02-22"

Apreciamos que tenemos un un dato diario entre el 26 de octubre y el 22 de febrero.

lapply(Daiyli_avgs3, function(x, titulo){
   k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
  ggplot(data=x, aes(x=Date))+
    geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
    geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+
    theme(panel.background = element_blank(),panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
                                           name = " LLuvia acumulada [mm]", 
                                           breaks = seq(min(cumsum(x$rDHI_acum)),
                                                        max(cumsum(x$rDHI_acum)),
                                                        by=50)),
                       breaks = seq(min(x$rDHI_acum),
                                    max(x$rDHI_acum),
                                    by=10))+
    ggtitle(titulo) +
    theme(plot.title = element_text(hjust = 0.5))

}, titulo=names(Daiyli_avgs3))
## $`-8.05770874023438__42.6729469299316`

## 
## $`-7.32501220703125__42.6962661743164`

## 
## $`-7.88290405273438__42.8138198852539`

## 
## $`-8.04904174804688__42.5383224487305`

Es interesante apreciar como pese a las altas correlaciones existe una pequeña diferencia al final del periodo en lo que a lluvia acumulada se refiere. Esto nos da pie a investigar cuales son los puntos que calculan mejor la lluvia acumulada.

vec_diff<- as.vector(unlist(lapply(Daiyli_avgs2, function(x){
  r<- sum(x$rWRF_acum)-sum(x$rDHI_acum)
  return(r)
})))

vec_diff[order(abs(vec_diff))]
##  [1]    4.404910    4.729176  -10.388001  -20.714521   31.905542
##  [6]  -32.582030   36.890532   39.813112  -49.308699  -49.663263
## [11]  -55.735768   76.225877  -77.398381  -79.743920  -81.969871
## [16]  -87.654145  -91.766186   91.816973  -98.120085  102.573777
## [21] -108.333345 -108.674158 -111.758937 -117.214062 -120.455026
## [26] -127.220242 -130.881619 -135.337830 -136.514876 -147.815890
## [31]  149.886050 -155.839446 -160.733827 -171.574046  180.299144
## [36] -185.876094  191.348111 -193.420657 -213.624562 -220.962254
## [41] -225.314246 -229.023156 -256.119342
Cor_rain_place_ccf[order(abs(vec_diff))[1:10],c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]

Vemos que entre las dos tablas que hemos sacado… la localización que se repite es la número 14, que se encuentra a 29 km. Esto es muy interesante para realizar luego las predicciones en base a esta localización.

Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10],c("AWRF_aport3","AWRF_ADHI", "Dist")]

Vemos que los puntos que se repiten en ambas listas… buscando mejor correlación con aportación, cercania y mejor resultado con la lluvia acumulada.

best_points<- order(abs(vec_diff))[1:10][order(abs(vec_diff))[1:10]%in%order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10]]

best_points
## [1] 28 18 19 32
litros_acum_diff<- vec_diff[best_points]
cbind(litros_acum_diff,
Cor_rain_place_ccf[best_points,c("AWRF_aport3","AWRF_ADHI", "Dist")])

Para un futuro tendremos que elegir entre estos 3 puntos. Añadimos también datos de correlación.

Daiyli_avgs4<-lapply(Daiyli_avgs2[c(14,18,19)], function(x)x[, c("Date","rWRF_acum",
                                      "aport_mean_SMA_diario", 
                                      "rDHI_acum")])
lapply(Daiyli_avgs4, function(x, titulo){
   k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
  ggplot(data=x, aes(x=ymd(Date)))+
    geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
    geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+
    theme(panel.background = element_blank(),panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
                                           name = " LLuvia acumulada [mm]", 
                                           breaks = seq(min(cumsum(x$rDHI_acum)),
                                                        max(cumsum(x$rDHI_acum)),
                                                        by=50)),
                       breaks = seq(min(x$rDHI_acum),
                                    max(x$rDHI_acum),
                                    by=10))+
    ggtitle(titulo) +
    theme(plot.title = element_text(hjust = 0.5))+
    geom_line(aes(y=aport_mean_SMA_diario/k), col="green",alpha= 0.9)

}, titulo=names(Daiyli_avgs4))
## $`-8.04904174804688__42.5383224487305`

## 
## $`-7.31787109375__42.5615844726562`

## 
## $`-8.05770874023438__42.6729469299316`

6 Predictive Model

Aquí viene lo gordo. Ya sabemos cuales son los puntos de WRF que mejor correlan con las observaciones aportadas por DHI. Lo primero, será cojer los datos horarios de los puntos anteriormente vistos. Cositas interesantes que quiero probar:

  • Se me acaba de ocurrir que igual se podría correlar guay la media movil de la precipitación con la aportación. Debido a esa inercia de la que hablamos

  • Tambien será necesario analizar como se comportar la variación del nivel con la aportación. en principio deberían correlar bien.

path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2018/10/01")),])
})


fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/12/01")
    
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})


cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
  y<- lapply(x, function(r){
    
    
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
}, 
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})

6.1 comprobación correlacion diferencia de nivel y aportacion

Realizamos una gráfica de la diferencia de nivel respecto a la aportación.

x<- cut_train[[1]]
plot(diff(x$nivel_mean_SMA)*20,  
     ylim = c(-1.5,2), 
     type = "l", 
     xlab = "",
     xaxt="n",
     yaxt="n",
     ylab = " ")
lines(x$aport_mean_SMA/100-1, col="red")

Vemos que existe un disparo en la gráfica… eso es porque existe un salto temporal de 3 días en diciembre… en los que el nivel vario 80 cm… cortamos los datos para eliminar ese disparo.

xx<- x[which(x$Date> ymd("2018/12/10")),]  
plot(diff(xx$nivel_mean_SMA)*20,  
     ylim = c(-1.5,2), 
     type = "l", 
     xlab = "",
     xaxt="n",
     yaxt="n",
     ylab = " ")
lines(xx$aport_mean_SMA/100-1, col="red")

Comprobamos las cross-correlation de estas dos variables.

ccf_nivel_aport<- ccf(diff(xx$nivel_mean_SMA),
                      xx$aport_mean_SMA,
                      lag.max = 1000)

print(paste0("Correlación máxima: ",
             max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
             ccf_nivel_aport$lag[
               which.max(ccf_nivel_aport$acf)], " horas"))
## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"

La aportación se ve reflejada en el nivel con un retraso de 15 horas con una correlación bastante alta de casi el 80%.

vec_cor<- vector()

for (i in seq(6, 144, by=2)) {
  x<- cut_train[[1]]
  x$prep_hourlySMA<- SMA(x$prep_hourly, i)
  x<- x[complete.cases(x),]

  rain_aport<- ccf(x$prep_hourlySMA, 
      x$aport_mean_SMA, 
      lag.max = 1000, 
      plot = F)
  
  vec_cor[i]<- max(rain_aport$acf)
  
  
}

plot(vec_cor, 
     xlab = "SMA points",
     ylab = "Corr")

Se ha decidido que el SMA ha emplear será de 48 puntos. A continuación realizamos una predicción de la aportación empleando diferentes métodos y variables.

  • Aportación de hace 3 días
  • Lluvia de hace 1 día
  • Lluvia actual
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2018/10/01")),])
})


fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/12/01")
    
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})


cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
  y<- lapply(x, function(r){
    
    
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
}, 
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})


x<- cut_train[[1]]
x$aport1<- lag(x$aport_mean_SMA, 24)
x$aport2<- lag(x$aport_mean_SMA, 48)
x$aport3<- lag(x$aport_mean_SMA, 72)
x$aport4<- lag(x$aport_mean_SMA, 96)
x$aport5<- lag(x$aport_mean_SMA, 120)

x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)

x<- x[complete.cases(x),]


y<- cut_predict[[1]]
y$aport1<- lag(y$aport_mean_SMA, 24)
y$aport2<- lag(y$aport_mean_SMA, 48)
y$aport3<- lag(y$aport_mean_SMA, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)

y<- y[complete.cases(y),]

indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")

y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]




fit_1  <- svm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2  <- svm(aport_mean_SMA  ~ prep_hourlySMA48_lag + aport3 +  prep_hourlySMA48, data = x)
fit_3  <- svm(aport_mean_SMA  ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)

fit_4  <- lm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5  <- lm(aport_mean_SMA  ~ prep_hourlySMA48_lag + aport3 + prep_hourlySMA48, data = x)
fit_6  <- lm(aport_mean_SMA  ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)


uncorrected<-y$aport_mean_SMA
pred_aport1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))

pred_aport4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))


plot(uncorrected, x=y$Date, type = "l",
     ylim = c(0,600), 
     xlab = paste0(range(y$Date)),
     ylab = "Aportación [m³/s]",
     main =  "Predicción de aportación \n empleando diferentes métodos")
lines(pred_aport1,x=y$Date, col="red")
lines(pred_aport2,x=y$Date, col="green")
lines(pred_aport3,x=y$Date,col="blue")
lines(pred_aport4,x=y$Date, col="red", lty=2)
lines(pred_aport5,x=y$Date, col="green", lty=2)
lines(pred_aport6,x=y$Date,col="blue", lty=2)

Del anterior gráfico podemos concluir que no vale la pena realizar la predicción aplicando svm (Support Vector Machine), o por lo menos no de esta forma…. Es interesante remarcar que solamente estamos en disposición de realizar un pronóstico a 3 días… dado que el modelo necesita el dato de aportación de hace 3 días para funcionar.

ccf(diff(x$nivel_mean_SMA), x$prep_hourlySMA48, lag.max = 1000)

print(paste0("Correlación máxima: ",
             max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
             ccf_nivel_aport$lag[
               which.max(ccf_nivel_aport$acf)], " horas"))
## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"

No tenemos valores de aportación en la página web, solamente nivel… por lo tanto hay que conseguir una buena correlación entre nivel y aportación. O directamente predecir el nivel con la lluvia.

path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r){ 
    r<- r[which(r$Date > ymd("2018/10/01")),]
  r$nivel_mean_SMA1<- c(diff(r$nivel_mean_SMA),0)
  return(r)
  })
})


fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/12/01")
    
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
  
  
})


cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
  y<- lapply(x, function(r){
    
    
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
}, 
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
  
  
})


x<- cut_train[[1]]
x$nivel1<- lag(x$nivel_mean_SMA1, 24)
x$nivel2<- lag(x$nivel_mean_SMA1, 48)
x$nivel3<- lag(x$nivel_mean_SMA1, 72)
x$nivel4<- lag(x$nivel_mean_SMA1, 96)
x$nivel5<- lag(x$nivel_mean_SMA1, 120)

x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)

x<- x[complete.cases(x),]


y<- cut_predict[[1]]
y$nivel1<- lag(y$nivel_mean_SMA1, 24)
y$nivel2<- lag(y$nivel_mean_SMA1, 48)
y$nivel3<- lag(y$nivel_mean_SMA1, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)

y<- y[complete.cases(y),]

indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")

y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]




fit_1  <- svm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2  <- svm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag + nivel3 +  prep_hourlySMA48, data = x)
fit_3  <- svm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)

fit_4  <- lm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5  <- lm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag + nivel3 + prep_hourlySMA48, data = x)
fit_6  <- lm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)


uncorrected<-y$nivel_mean_SMA1
pred_nivel1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))

pred_nivel4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))


plot(uncorrected, x=y$Date, type = "l",ylim=c(0, 0.2),
     xlab = paste0(range(y$Date)),
     ylab = "nivelación [m³/s]",
     main =  "Predicción de nivelación \n empleando diferentes métodos")
lines(pred_nivel1,x=y$Date, col="red")
lines(pred_nivel2,x=y$Date, col="green")
lines(pred_nivel3,x=y$Date,col="blue")
lines(pred_nivel4,x=y$Date, col="red", lty=2)
lines(pred_nivel5,x=y$Date, col="green", lty=2)
lines(pred_nivel6,x=y$Date,col="blue", lty=2)